Isotropic singularity in inhomogeneous brane cosmological models 
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We discuss the asymptotic dynamical evolution of spatially inhomogeneous brane-world cosmo- 
logical models close to the initial singularity. By introducing suitable scale-invariant dependent 
variables and a suitable gauge, we write the evolution equations of the spatially inhomogeneous Gi 
brane cosmological models with one spatial degree of freedom as a system of autonomous first-order 
partial differential equations. We study the system numerically, and we find that there always exists 
an initial singularity, which is characterized by the fact that spatial derivatives are dynamically neg- 
ligible. More importantly, from the numerical analysis we conclude that there is an initial isotropic 
singularity in all of these spatially inhomogeneous brane cosmologies for a range of parameter values 
which include the physically important cases of radiation and a scalar field source. The numerical 
results are supported by a qualitative dynamical analysis and a calculation of the past asymptotic 
decay rates. Although the analysis is local in nature, the numerics indicates that the singularity is 
isotropic for all relevant initial conditions. Therefore this analysis, and a preliminary investigation 
of general inhomogeneous (Go) models, indicates that it is plausible that the initial singularity is 
isotropic in spatially inhomogeneous brane-world cosmological models and consequently that brane 
cosmology naturally gives rise to a set of initial data that provide the conditions for inflation to 
subsequently take place. 

I. INTRODUCTION 

String inspired theories in which the matter fields are confined to a 3-dimensional 'brane-world' embedded in 
1 + 3 + d dimensions while the gravitational field can also propagate in the d extra dimensions [1] are currently of 
great interest. In particular, Randall and Sundrum [2] have shown that for d = 1, gravity can be localized on a 
single 3-brane at lower energies even when the fifth dimension is infinite. The Randall-Sundrum type models are 
simple phenomenological models, which capture some of the essential features of the dimensional reduction of eleven- 
dimensional supergravity introduced by Hof ava and Wittcn [3] . An elegant geometric formulation and generalization 
of the Randall-Sundrum-type brane-world models has been given [4,5]. Recently there has been great interest in 
Randall-Sundrum-type brane-world cosmological models [2], particularly in an attempt to understand the dynamics 
of the universe at early times. Brane-world models have a different qualitative behaviour than their general-relativistic 
counterpart [6,4], especially at high energies when the energy density of the matter is larger than the brane tension 
and the behaviour deviates significantly from the classical case. 

The asymptotic dynamical evolution of spatially homogeneous brane-world cosmological models close to the initial 
singularity was studied in [7]. It was found that an isotropic singularity [8] is a past-attractor in all orthogonal 
Bianchi models and is a local past-attractor in a class of inhomogeneous brane-world models (and consequently these 
models do not exhibit Mixmaster or chaotic- like behaviour close to the initial singularity). However, the study of 
the behaviour of spatially homogeneous brane-worlds close to the initial singularity in the presence of both local and 
nonlocal stresses indicates that for physically relevent values of the equation of state parameter there exist two local 
past attractors for these brane-worlds, one isotropic past attractor and one anisotropic past attractor; i.e., in these 
brane-worlds the initial singularity can be locally either isotropic or anisotropic [9,10] (however, the anisotropic models 
appear to be unphysical and can likely be ruled out). Therefore, it is plausible that typically the initial singularity 
is isotropic in the brane-world scenario. Consequently, it was suggested that brane cosmology naturally gives rise to 
a set of initial data that provide the conditions for inflation to subsequently take place, thereby solving the initial 
conditions problem and leading to a self-consistent and viable cosmology [11]. We argued [7] that it is plausible that 
typically the cosmological singularity is isotropic in spatially inhomogeneous models. We shall study this further here. 

We shall consider the dynamics of a class of spatially inhomogeneous cosmological models with one spatial degree of 
freedom in the brane-world scenario. The G2 cosmological models admit a 2-parameter Abelian isometry group acting 
transitively on spacelike 2-surfaces. These models admit one degree of freedom as regards spatial inhomogeneity, and 
the resulting governing system of evolution equations constitute a system of autonomous partial differential equations 
in two independent variables. We follow the formalism of [12] which utilizes area expansion normalized scale-invariant 
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dependent variables, and we use the timclike area gauge to discuss the asymptotic evolution of the class of orthogonally 
transitive G2 cosmologies near the cosmological initial singularity. In this article we shall consider numerically the 
local dynamical behaviour of this class of spatially inhomogeneous models close to the singularity. 



A. Governing equations 

The field equations induced on the brane, using the Gauss-Codazzi equations, matching conditions and Z2 symmetry, 
result in a modification of the standard Einstein equations with the new terms carrying bulk effects onto the brane 
[4,5]: 

Gfx V = K T^ v + K Sfj^u — £fj_ u , (1) 

where 

k 2 = 8tt/M p 2 , A = 6^J (2) 

For any matter fields (scalar fields, perfect fluids, kinetic gases, dissipative fluids, etc.), including a combination of 
different fields, the general form of the brane energy-momentum tensor can be covariantly given as 

= pu^Ut, + ph^ + ti> + 2q {}1 u v) . (3) 

The decomposition is irreducible for any chosen 4- velocity . Here p and p are the energy density and isotropic 
pressure, and h^ v = g^ v + u^u v projects orthogonal to u^. The energy flux obeys = q^), and the anisotropic stress 
obeys ir^ — ^(^ v ), where angle brackets denote the projected, symmetric and tracefree part: 

v<„) - K v v v , w {llv) = [h (ll a h v f - \h af) v] w aP ■ 

We shall primaily be interested in perfect fluid sources obeying the linear equation of state p = (7 — 1) p, where 7 is 
constant with < 7 < 2. The case 7 = 4/3 corresponds to radiation. Scalar fields correspond to a stiff fluid 7 = 2 
close to the initial singularity. 

The dynamical equations on the 3-brane differ from the general relativity equations [4,5] in that there are nonlocal 
effects from the free gravitational field in the bulk, transmitted via the projection £ M „ of the bulk Weyl tensor, and 
the local energy-momentum corrections, which are significant at very high energies and particularly close to the initial 
singularity. The matter fields contribute local quadratic energy- momentum corrections via the tensor S^ u , given by 

Sfiv "yI^ol — -jT^ a T „ + 2i9fif 3T a pT ^ {T a ) . (4) 

Equations (3) and (4) imply the irreducible decomposition of 1 

= ± [2p 2 - 3ir aP Tr afj ] u^u, + £ [ 2/9 2 + 4pp + ir a pir af) - Aq a q a ] 

- j2-(p+ Sp)TT^ - \-K a ^TX v) a + \q^q v ) + \pq{ p u v) - \q a ir a{p ,u v) . (5) 

The quadratic energy-momentum corrections to standard general relativity will be significant for K 4 p 2 > n 2 p in the 
high-energy regime close to the singularity. 

Nonlocal effects from the bulk can be irreducibly decomposed into 

£iiv = - (J^J [U ( V + \hy, v ) + + 2Q( li u v - ) ] , (6) 

in terms of an effective nonlocal energy density on the brane, U, arising from the free gravitational field in the bulk, 
an effective nonlocal anisotropic stress on the brane, V^, and an effective nonlocal energy flux on the brane, <2 M [5]. 



1 We take this opportunity to correct equation (7) of [5]. 
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All of the bulk corrections may be consolidated into an effective total energy density, pressure, anisotropic stress 
and energy flux, as follows. The modified Einstein equations take the standard Einstein form with a redefined 
energy-momentum tensor: 



where 
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(Note that k 4 /k 6 is dimensionless.) 

As a consequence of the form of the bulk energy-momentum tensor and of Zi symmetry, it follows [4] that the 
brane energy-momentum tensor separately satisfies the conservation equations, i.e., 

VT,„=0. (14) 

Consequently, the Bianchi identities on the brane imply that the projected Weyl tensor obeys the non-local constraint 

WE^ = k 4 V»S^ . (15) 

The results of [7] are incomplete in that a description of the gravitational field in the bulk is not provided. Unfor- 
tunately, the evolution of the anisotropic stress part is not determined on the brane. The correction terms must be 
consistently derived from the higher-dimensional equations. Since V^v corresponds to gravitational waves in higher- 
dimensions, it is expected that the dynamics will not be affected significantly at early times close to the singularity 
(see [13] and later). Henceforth we shall effectively assume that 



(16) 



When V^u = 0, the evolution of £ M „ is fully determined [4]. In the inhomogeneous cosmological models of interest 
here a non-zero D^p acts as a source for Q M , and hence = is not consistent with an inhomogeneous energy 
density, and we need to include a dynamical analysis of the evolution of Q M . We shall make no further assumptions 
on the models and include all terms in the numerical analysis. 



B. Initial Singularity 

From the numerical analysis we shall find that the area expansion rate increases without bound (f3 — ► oo) and the 
normalized frame variable [12] vanishes (Ei 1 — > 0) as logarithmic time t — ► — oo. Since [3 — > oo (and hence the Hubble 
rate diverges), there always exists an initial singularity as t — > — oo. Thus the singularity is characterized by Ei 1 = 0, 
which allows both dynamical and numerical results to be obtained (see later). 

In [7] it was shown that the total energy density p^ooast^— oo. It then follows directly from the conservation 
laws that pb ~ p\ dominates as t ^ — oo and that all of the other contributions to the brane energy density are 
negligible dynamically as the singularity is approached. The fact that the effective equation of state at high densities 
become ultra stiff, so that the matter can dominate the shear dynamically, is a unique feature of brane cosmology. 
Hence close to the singularity the matter contribution is effectively given by 
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2A' 

-tat = 1 (p2 + 2 pp) = (2 7 - 1)A 6 . (18) 



In addition, it follows from the conservation laws that f2 — > 0, fi^ — > and Oa ^ as the initial singularity is 
approached. Hence, the models isotropize to the past. We shall study the generality of this result. 

Models with an isotropic initial singularity [8] satisfy lim t ^_ 00 Ct b = 1, lim t ^_ oc u = 0, lim t ^_ 00 S 2 = 0. Their 
evolution near the cosmological initial singularity is approximated by the flat model corresponding to the 'equilibrium 
point' Th-, characterized by 2 

n b = i; o = ^i 1 = £+ = s_ = s x = n_ = a x = v = q u = n = n u . (19) 

Tb corresponds to a spatially homogeneous and isotropic non-general-relativistic brane-world model (which is valid 
at very high energies (p ^> A) as the initial singularity is approached). Note that these solutions are self-similar, and 
are referred to as Binctruy, Dcffayet and Langlois solutions [6] or Brane- Robertson- Walker models [11]. It was shown 
that for all physically relevant values of 7 (7 > 1), T b is a source (or past-attractor), and hence the singularity is 
isotropic, in non-tilting spatially homogeneous brane-world models [7]. It was also shown that T b is a local source or 
past-attractor in the family of spatially inhomogeneous 'non-tilting' G2 cosmological models for 7 > 1 [7] . 

In this paper we shall study the nature of the initial singularity in spatially inhomogeneous brane cosmological 
models. In particular, we shall study numerically the class of G2 models. An analysis of the behaviour of spatially 
inhomogeneous solutions to Einstein's equations near an initial singularity has been made in classical general relativity; 
in an investigation of a class of Abelian G2 spatially inhomogeneous models [14] and a numerical investigation of a 
class of vacuum Gowdy G2 cosmological spacetimes [15], it was shown that the presence of the inhomogeneity ceases 
to govern the dynamics asymptotically toward the singularity. 

II. BRANE G 2 COSMOLOGY 

We shall consider the class of G2 cosmologies with two commuting Killing vector fields, which consequenly admit 
one degree of spatial freedom [16]. We shall follow the approach of van Elst, Uggla and Wainwright [12]. The evolution 
system of the EFE are partial differential equations (PDE) in two independent variables. The orthonormal frame 
formalism is utilized [17,18] with the result that (i) the governing equations are a first-order autonomous equation 
system, (ii) the dependent variables are scale- invariant. In particular, we define scale-invariant dependent variables 
by normalisation with the area expansion rate of the G2 _ orbits in order to obtain the evolution equations as a system 
of PDE, ensuring the local existence, uniqueness and stability of solutions to the Cauchy initial value problem for G2 
cosmologies. Following [12] we assume that the Abelian G2 isometry group acts orthogonally transitively on spacelike 
2-surfaces, and introduce a group- invariant orthonormal frame {e Q }, with e2 and e3 tangent to the C?2-orbits. The 
frame vector field eo, which defines a future- directed timelike reference congruence, is orthogonal to the G2-orbits and 
it is hypersurface orthogonal and hence is orthogonal to a locally defined family of spacelike 3-surfaces t = const. We 
then introduce a set of symmetry-adapted local coordinates { t, x, y, z } 

e = N- 1 d t , e 1 = e 1 1 8 x , e A = e A B d x s , A,B = 2,3, (20) 

where the coefficients are functions of the independent variables t and x only. The only non-zero frame variables are 
thus given by 

N , ex 1 , e A B , (21) 

which yield the following non-zero connection variables [18]: 

a, (3, 01, n+, a-, n x , cr x , n_, u\, Vt\. (22) 

The variables a, (3, cr_ and cr x are related to the Hubble volume expansion rate H and the shear rate o~ a p of the 
timelike reference congruence e ; in particular, 6 := 3H = a + 2/3. The variables a\, n+, n x and n_ describe the 



2 A11 of the variables used here are defined later (e.g., see equations (67)-(69)). 



non-zero components of the purely spatial commutation functions a a and n a p [16]. Finally, the variable u\ is the 
acceleration of the timelike reference congruence eo, while 17 1 represents the rotational freedom of the spatial frame 
{e a } in the (e 2 , e 3 )-plane. Setting fli to zero corresponds to the choice of a Fermi-propagated orthonormal frame 
{ e a }. Within the present framework the dependent variables 

{W.ui.fii} (23) 

enter the evolution system as freely prescribable gauge source functions. 

Since the G2 isometry group acts orthogonally transitively the 4-velocity vector field u of the perfect fluid is 
orthogonal to the G2-orbits, and hence has the form 

u = r(e +we 1 ), (24) 

where T= (1 - v 2 )" 1 / 2 . 

We assume that the ordinary matter is a perfect fluid with equation of state 

Pfl = (7 - ±)pfl • (25) 

In a tilted frame, we have 

T^u = pu^u v + ph^ v + it + 2q^u v ) , (26) 

where 

G+ 1 3(7-l)(l-« 2 )+7" 2 IP IP (07 , 

P = Y~^ Pfu p = 3 Gl p ' qa = G^ Va, na0 = G^ v ( aV P) { - 27 > 

and G+ = 1 + (7 — l)v 2 . The basic variables that we use are p and v a . 
The quadratic correction matter tensor S a p is given by 



where 



•V = P bu ^ + P b V + *V + 2( l\n u v) > (28) 



p* = -L(2p 2 -3n a ^) (29) 



1 (l-v 2 ) 
12 G+ 2 



[1 + (2 7 - l)v 2 ]p 2 , (30) 



p b = ^(2p 2 + 4pp + 7T aP 7T^ -4q a q a ) (31) 
_ 13(2 7 -l)(l- w 2) + 2 7 ^^ 



3 1 + (27 - l)v 
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q b a = ^(2pq a -3q^ a ) (33) 



2 7 p" 



7T b 



(34) 



1 + (2 7 - l)v 2 a ' 

— [-(p + 3p)7r Q/3 - 37r 7(Q 7r (3) 7 + 3^^] (35) 
2 7 p fc 

U(a^/3> • (36) 



1 + (27 - l)w 2 

Comparing with (27), we see that the quadratic matter is effectively a perfect fluid with equation of state 

b 

fl ■ 

From equation (6), the bulk matter E^y is given by 



P b fl = (2 7 -l)p b fl . (37) 



6k 2 

£^ = - — (fi u u^u v +p u h^ + ttV + 2q u {fl u v) ) (38) 
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where 



1 1 + \v 2 

p u = ^ = T^P u n (39) 



^ = ^=iT^^ + 7^ {qUf ^ (41) 
i -p u 

n%„ = Vi V, v = ^j^v^ + (ttW . (42) 
Comparing with (27), we see that the bulk matter is similar to a perfect fluid with equation of state 

P u fi = \p u fi , (43) 

but with non-zero (q u fi)^ and (w u fi)p V - We shall use p u and q"^ as the basic variables, and set (tt u ji)^ — in this 
paper (see (16)). 

Since 

9% = 9"<^\ Vf, = v5^, (tt 11 /;)^ = 0, 
we can define q, qb, q u , tt, ttj, and ir u as follows: 

(91, 92, 9 3 ) = (9, 0, 0), (<A, q b 2 , q b 3 ) = (q b , 0, 0), {q\, q u 2 , q u 3 ) = (q u , 0, 0) 

and 

tt^u = when /i 7^ v. ttu = n, 7^2=^33 = -^ = -^ 

6 n V, -L b b b 7I ' b ll n b 

7T ^ = when ^ T V- n°u = n b , tt 22 = ?r 33 = — = — — 

= when /j/k tt n = tt„, tt 22 = tt 33 = — = — — 

The orthonormal frame version of the EFE, matter equations and non-local equations, when specialised to the 
orthogonally transitive Abelian G2 case with the dependent variables presented above [12], takes the following form: 

Einstein field equations and Jacobi identities 

Evolution equations: 

N^ 1 d t a = - a 2 + [3 2 - 3 (a 2 _ - n\ + a\ - n 2 _) -a\ + (e^ d x + in) u x 

-i(p tot +p tot )+7r tot (44) 
iV" 1 d t p = - f /3 2 - § (a 2 _ + n\ + a 2 , + n 2 _) - \ (2«i - 01) 01 

-i(p tot +^ tot ) (45) 

TV" 1 d tai =-(3{ii 1 + a 1 )-3 (n x <7_ - n_ a x ) - \q tot (46) 

TV -1 dtn + = — an + + 6 (<r_ n_ + a x n x ) — (ei 1 d x + in) f^i (47) 

A^ 1 dta- + ex 1 d x n x = - (a + 2/3) a- - 2 n+ n_ - (in - 2a x ) n x - 2 Oi cr x (48) 

TV -1 <9 t n x + ei 1 <9 X <7_ = — an x + 2 a x n + — ii\ cr_ + 2 Oi ra_ (49) 

A^ 1 cVx - ei 1 9 x n_ = - (a + 2/3) a x - 2 n+ n x + (in - 2oi) n_ + 2 Q x cr_ (50) 

/V -1 <9 t n_ — ei 1 <9 X <7 X = — a + 2 cr_ ra + + ii\ a x — 2 Oi n x . (51) 

Constraint equations: 
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0= 2(2e 1 1 d x -3ai)ai - 6 (n 2 + n 2 _) + 2 (2a + /?) /3 - 6 (cr 2 + _ 2p tot (52) 
0= ei 1 9,/? + oi (a - /?) - 3 (n x a_ - n_ (j x ) - ^ tot . (53) 

Bianchi identities (conservation equations) 

Evolution equations: 

— i N ~' dt + ^v ei 1 d x ) p + h ei 1 d x v = -^- /i [ a (1 + w 2 ) + 2/3 + 2 (tii - at) v ] (54) 

p G+ G + 

^ p (iV- 1 9, - - A_ « ei' d *) v + h ei 1 d xP = - p(l-v 2 )l(2- 1 )av-2( 1 -l)(3v 

Jl G + G_ JlCr- 

+ G_ wi + 2 ( 7 - 1) ai w 2 ] , (55) 



where 



h ■= ^ (1 - v 2 ) 2 , h:=^(l-v 2 ) 2 , (56) 
G± ~l±( 1 -l)v 2 , / 3 :=(3 7 -4)-( 7 -l)(4- 7 ) W 2 . (57) 



From equation (15) we obtain: 
Non-local conservation equations 

Evolution equation: 



N- 1 d t pu + e l 1 d x q u = ~ ~ ^ [(1 + v 2 )a + 2(3} p u - 2(ui - oi)g„ + (58) 

3 + ?r 

A^ 1 a t g u + ei 1 a x (p u + 7T„) = [(1 + « 2 )"i - 2« 2 ai] p u - 2(a + (3)q u + Y (59) 

3 + v z 



where 



l + 3v 2 16v 
3 + ir (3 + ?r) z 

7 2 w(l — v 2 )p 2 



d x (p u + ttu) = — — ^d x p u + 2 , 2 p u d x v (60) 



Y = ' 6 a + 2 G [(l-« a )a + 2fl-2H 

7(1 - ^ 2 ) V i g fl | 7 2 ^(1 - ^ 2 )(3 + (7 - IK )P 2 ig , . 

" 6G + 2 G _ 61 + 6G+^ 61 ^ ' (61) 



and 



tot 2 6k2 6k 2 

J tot = K 2 p + — Pb + — Pv 



= * p + ^l2G^ i 1 + ( 2 "'- l > 2 ) P 2 + Jfc ( 62 ) 

p t0t = K 2 p + — Pb + — Pu 

1 3(7-l)(l-^ 2 )+7^ 2 , r2l , 6 K 2 1 3(2 7 -l)(l- t ; 2 ) + 27i; 2 2 2 6« 2 1 

= 3 G+ KP+ — 3 ( )P + ~3^ (63) 

tot 2 6k2 6k2 

9 = K 2 q + — q b + —q u 



^^^^ 



tot 2 , 6k 2 6k 2 

7T = K 7T H — TT b H — TT V 

A A 



2 7 t; 2 2 6re 2 2-fv 2 (1 - w 2 ) , 6re 2 2 4w 2 

3G7 K ^+~3 6G + 2 P + — 3(37^"° (65) 
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Finally, we have: 
Gauge fixing condition: 

= N- 1 e 1 1 d x N-ii 1 . (66) 
The frame variables eA B decouple from the remaining equations and we need not consider them further. 

A. Scale-invariant reduced equation system 

We introduce /3-normalised frame, connection and curvature variables as follows [12] 3 : 

(TV -1 , Ei 1 ) 

(U, A, (1-3S+), E_, N x , S x , N_,N+, R) 

( ^1 1 Qu ) 

B 

where (3 is the area expansion rate of the G 2 -orbits. The new dimensionless dependent variables are invariant under 
arbitrary scale transformations, and are linked to the _ff-normaliscd variables through the relation H = (1 — E + ) (3 
[7,16]. Note that in the units we have chosen the matter variable v is already dimensionless. 

In order to write the dimensional equation system in terms of scale-invariant dependent variables (67) - (69) it 
is necessary to introduce the time and space rates of change of the normalisation factor /3. We use the evolution 
equation (45) and the constraint equation (53) to define q and r in terms of the remaining scale-invariant dependent 
variables: 

N- 1 d t B = -(q+l)B (71) 
0=(E 1 1 d x + r)B , (72) 

the expressions (90) and (91) (below) for q and r are purely algebraical, and are referred to as the defining equations 
for q and r; q and r play the role of an "area deceleration parameter" and a role analogous to a "Hubble spatial 
gradient". Using equations (71) and (72), the definitions (67) - (69) and equation (13) in [12], it is straightforward 



to transform the dimensional equation system to a /3-normalised dimensionless form. 
Scale-invariant equation system 

Evolution system: 

M- 1 d t B = -{q+l)B (73) 

TV -1 dtEx 1 = (q + 3S + ) Ei 1 (74) 

AT 1 d t ^+ = (q + 3S+ - 2) S+ - 2 (N 2 + N 2 X ) - ± E, 1 d x r - § II tota] (75) 

N- 1 d t A= (<? + 3E+) A + (r - {/) (76) 

AT 1 d t N+ = {q + 3E+) N+ + 6 (£_ N. + S x N x ) - (V d x - r + U) R (77) 

Af- 1 d t Z- + E! 1 d x N x = (q + 3S+ - 2) £_ - 2 N+ 7V_ + (r - U + 2A) N x -2i?FJ x (78) 

JV _1 $JV X + E 1 1 d x Y,_ = (q + 3Z+)N x + 2E XJ /V + + (r-f/)£_ + 2RN_ (79) 

M~ 1 d t T, x -E 1 l d x N_ = (q + 3S+ -2)E X - 2 N+ N x - (r - U + 2A) N_ + 2 i?FJ_ (80) 

TV -1 d t N_ — Ei 1 d x Yi x = (q + 3E+) 7V_ + 2 £_ 7V+ - (r - L>)S X -2RN X . (81) 



§ (AT 1 9 t + -f- w^ 1 8 X ) fi + h Ei 1 d x v = 2 1£ [ -+ (q + 1) - \ (1 - 3E+) (1 + v 2 ) 



3 Expressions for the area expansion rate and the scale-invariant dependent variables in terms of the line element, written in 
the separable area gauge, are given in [12]. 



(TV-Vr 1 )//? 

(iti, ai, a, n x , cr x , n_, rc+, Qi )//3 

— P«, — (MAW = ^S-p, Afi . 2 



(n 2 p, ^- Pu , ^q u )/(3(3 2 ) = (K 2 p, — , -4)/(3/3 2 ) 



(67) 
(68) 

(69) 
(70) 
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— l + (r — U + A)v] 



^ n (Af- 1 dt - 7 ^ r - vE 1 1 d x ) v + h Ei 1 d x n - 2 -A- n (1 - v 2 ) [ ^LJl 

Ji G + G_ 7iCx_ 7 



(1 -v 2 )r 



1 + v 2 
2(,-l) + 12^E + 



AT 1 a t Q u + E^d^Pu + n„) = 2(g + 3S+ - 1)Q„ 



3 + v 



- |(2-7) (l-3E+)« 

+ (7-l)(l-.4t;)w- |G_f7 

ft„ - 2({/ - r - A)Q„ + 

2(1 + -y 2 )(f7 - r) - 4u 2 A + (1 - v 2 )r 



where 



1 + 3u 2 16u 

d x (P u + I1 M ) = — d x fl u + 77T- 5^7^u d x v 



X = 



3 + v 2 

^(l-v^B 2 ^ 2 



(3 + v 2 y 



6G+ 2 G_ 



7 w(l - v 2 )(l - 3£+) + 2 7 -y - 2jv 2 A + 2(1 - w 2 )r 

B 2 n ESd x n + 7Ml ~ " 2)(3 + (7 ~ 1)t,2) bW E^ dx v . 



Constraint equations: 



6G+ 2 G_ 6G+ 3 G_ 
£V <9 X A = 2S+ - 1 + E 2 . + £ 2 + iV 2 + iV 2 + | r A + A 2 + Q total 



2 p 1 

3 



= r+ 3AS+ + 3(iV x JV_£ X ) + | Q total 



Defining equations for q and r 

+ 

E^d x B 



q := i + i ( 2 {7 - A) A + § (E 2 _ + iV 2 + S 2 + iV 2 ) + §( P total + n total ) 



r := 



B 



where 



Ntotal — + fife + Q u 
(l-v 2 ) 
~ + 12G+ 2 

Ptotal = P + Pb + Pu 



l + (2 7 -l)w 2 ) Q 2 B 2 +Q U 



3(f-l)(l-v 2 )+jv 2 3(2 7- 1)(1 - v 2 ) + 2 7 -, 



3G+ 

Qtotal = Q + Qb + Qu 

~G+ + 6G+ 2 
n to tai = n + n 6 + n u 



-0 + 



36 G+ 5 



-(l- v 2 )ft 2 B 2 + ifl 



ft 2 s 2 + 



2 7 ^ 2 7^ 2 (l-^) o2 p2 , 8i; 2 
3G+ + 9G+ 2 + 3(3 + « 2 ) 



define the various physical quantities. In particular, we note that 



n b ^^i(i + (2 1 -i)v 2 )n 2 B 2 . 



Gauge fixing condition: 



= AT 1 d x Af + (r - £7) . 



(82) 

(83) 

(84) 
(85) 

(86) 

(87) 

(88) 
(89) 

(90) 
(91) 



(92) 
(93) 
(94) 
(95) 

(96) 
(97) 



In the above, /i, / 2 , 7*3 and G± are defined by equations. (56) and (57), respectively. 
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B. Gauge choice 



The scale- invariant equation system in subsection II A contains evolution equations for the dependent variables 

{E, 1 , £+, A, N+, E_, N x , E x , JV_, ft 6 , «, Q} , (98) 
but not for the gauge source functions 

W t>, R} , (99) 

which are arbitrarily prescribable real-valued functions of the independent variables t and x, and thus does not 
uniquely determine the evolution of the G2 cosmologies. The reason for this deficiency is that the orthonormal frame 
{ e a } and the local coordinates { t, x } have not been specified uniquely. The remaining gauge freedom consists of a 
choice timelike reference congruence eo and of local time and space coordinates t and x (the temporal gauge freedom), 
and a choice of spatial frame vector fields e 2 and e3 (the spatial gauge freedom). 

We shall fix the spatial gauge by requiring 

N + = VSN_, i?=-V3E x , (100) 

which is preserved under evolution and under a boost. With this choice the evolution equation (77) becomes identical 
to equation. (81), and thus can be omitted from the full scale-invariant equation system. 

We fix the temporal gauge by adapting the evolution of the gauge source function U. The separable area gauge is 
determined by imposing the condition 

= (r - U) , (101) 

which determines U algebraically through equation. (91). There is thus no need to determine an evolution equation 
for U. It follows immediately from the gauge fixing condition (97) that Af = f(t). We now use the i-reparametrisation 
to set f(t) — No, a constant, which we choose to be unity, i.e., 

Af = AT :=l. (102) 

Therefore, t is effectively a logarithmic proper time, and the initial singularity occurs for t — > —00. 

The area density A of the G2-orbits plays a prominent role for G2 cosmologies. Expressed in terms of the coordinate 
components of the frame vector fields tangent to the G2-orbits this becomes 

A' 1 = e 2 2 e 3 3 - e 2 3 e 3 2 . (103) 

In terms of our scale-invariant dependent variables, the area density A of the G2~orbits satisfies the relations [12] 

A' 1 N- 1 d t A = 2 , A- 1 E 1 1 d x A = -2A . (104) 

Combining the two, the magnitude of the spacetime gradient V ' a A is 

(V a ^4) (V a ^4) = - 4/3 2 (1 - A 2 ) A 2 , (105) 

so V a *4 is timelike for A 2 < 1. 

For the class of G2 cosmologies in which the spacetime gradient V a A is timelike, we can choose the gauge condition 

A = , (106) 

which would be achieved by choosing e to be parallel to V a A It follows from equation. (104) that d x A = 0, and we 
obtain A = 1$ e 2J ^ at , which is function of t only (the so-called area time coordinate). This is rcfcrcd to as the timelike 
area gauge [12]. 

There are other gauge choices, such as the fluid comoving gauge or the synchronons gauge, but these are less 
convenient for numerical analysis. 



10 



C. Governing equations in timelike area gauge 



Let us explicitly give the evolution system in the timelike area gauge using equations (100), (101), (102) and (106). 
In the timelike area gauge we can use equation (100) to eliminate the evolution equation (77), and equation (76) 
becomes trivial. The relevant equations are: 
Evolution system: 



d t B = -(q + l)B 
dtE^ = ( <Z + 3S+)£; 1 1 
<9 t E + = (q + 3E+ - 2) E+ - 2 (Nl 



Nt 



■E^dxr- |n tota] 



d t E_ + E 1 1 d x N x 
d t N x +E 1 1 x ^ 
dtZx-E^dxN- 
dtN- - Eh 1 d x E x 



(q + 3E+ - 2) E_ + 2V3E 2 - 2V3 N 2 _ 
(<Z + 3E + )7V X 

(q + 3E+ - 2 - 2V3£_) E x - 2V3iV x N. 
(q + 3E+ + 2\/3E_) 7V_ + 2^3 E x N x 



^(d t + ^-vE 1 1 d x )Q + f 2 E 1 1 d x v = 2^-f 1 [^(q + l)- i(l-3£ + )(l+« 2 )- 1] 



f i n{dt ~ — 



vE 1 1 d x )v+ f 2 E 1 1 o x n = - 



h 



^(i-, 2 )[ (2 - 7) 



G+r 



fiG- 

+ (2- 7 )(l-3E + )t;-2(7-l)t;] 



1 + v 2 
2(q - 1) + 12^— -E+ 



3 + v 2 



n u + vx 



d t Qu + E^d^Pu + n„) = 2(q + 3E+ - 1)Q„ 



2(1 ~^ 2 ) 
3 + v 2 



where 



1 + 3u 2 16u 



Constraint equations: 



3 + v 2 
7 (l-v 2 )P 2 2 



2\2 



fl u d x v 



6G+ 2 G_ 
7 (1 -v 2 ) 2 
6G+ 2 G_ 



jv(l - v 2 )(l - 3E+) + 2jv + 2(1 - v 2 )r 



B 2 nE 1 1 d x fl + 



j 2 v(l ~ v 2 ){3 + (~f - l)v 2 ) 
6G+ 3 G_ 



P 2 f7 2 E^d^ 



= 2E+ - 1 + E 2 _ + E 2 + iV 2 + iV 2 + n + ^ (1 + (2 7 - 1) « 2 ) ft 2 P 2 



= r + 3(W x E_ -7V_E X ) + 



-n 



6G+ 2 



ft 2 B + 0, 



Defining equations for q and r: 

q := i + | (E 2 + iV 2 + E 2 + TV 2 ; 



27 - 1 + v2 (i - « 2 ) tt 2 b 2 + + 1t4±^ 



r := 



2 + 3(1- w 2 ) 

E X X 3 X B 
B 



12G+ 2 v 7 3 + v 2 

I" Ptotal + n to t£ 

( 7 -i)(i- w 2 )r! 2 i? 2 o u 



q + 3 E + — 2 + |( — fitotal + Ptotal + Iltotal) 



12 G+ 2 



3 + v 2 



G+ 

(2-7)" 
2G+ 



(107) 
(108) 
(109) 

(110) 
(111) 
(112) 
(113) 

(114) 



(115) 
(116) 

(117) 
(118) 

(119) 

(120) 
(121) 



(122) 



(123) 



(124) 



There is numerical (section III) and dynamical (section V) evidence for the existance of a number of monotonic 
functions close to the initial singularity. In particular, from equations (107) and (124) (using (37 — 1) > 0), close to 
an isotropic singularity [3 is itself monotonic. 
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III. NUMERICAL RESULTS 



We have written the governing equations as a system of evolution equations subject to the constraint equations 
(120) and (121). We can use (120) to obtain O u and (121) to solve for Q u and thus treat the governing system as 
a system of evolution equations without constraints (i.e., we don't need to use the evolution equations for Q, u and 
Q u in the numerics). In the numerical analysis we use the standard CLAWPACK package for PDEs with one space 
variable (see [19] for background). In the numerical calculations we prescribe periodic boundary conditions (we also 
implement Roe-averaging for the Riemann solver and choose Godunov splitting for source term splitting). 

From the numerical analysis we find that the area expansion rate increases without bound ((3 —* oo) and the 
normalized frame variable [12] vanishes (-Ei 1 — > 0) as t — > — oo. Since j3 — > oo (and hence the Hubble rate diverges), 
there always exists an initial singularity. In addition, we find that {ft, S_, N x , S x , N_, r} — > as t — > — oo for all 
7 > 1. In the case 7 > 4/3, the numerics indicate that {v, Q u , Qu, — > (and — * 1) for all initial conditions. 
In the case of radiation (7 = 4/3), the models still isotropize as t — > —00, albeit slowly. We shall investigate this 
degenerate case in more detail. For 7 < 4/3, {v, Q Ul £+} tend to constant but non-zero values as t — > —00. It is 
interesting to note that v 2 -/* 1; i.e., the tilt does not tend to an extreme value. 

The numerical results support the fact that all cosmological models have an initial singularity and that for the range 
of values of the equation of state parameter 7 > 4/3 the singularity is isotropic. Indeed, the singularity is isotropic for 
all initial conditions (and not just for models close to Tb) indicating that for 7 > 4/3 this is the global behaviour. We 
illustrate the numerical results for 7 = 1.8 by presenting some graphs of a typical run (see FIGs.1-3; in the TABLEs 
dlr denotes d x r) which show snapshots of the initial conditions and at earlier times, indicating isotropization. The 
numerical results support the exponential decay (to the past) of the anisotropics of section IV and [7] (see FIGs.4-5). 
Indeed, we found no evidence that models do not isotropize to the past for 7 > 4/3 [20]. 

We noted in [7], that an inhomogeneous energy density with a non-zero D u p acts as a source for Q^, and we must 
check that the evolution of Q M is consistent with the approximations used and consequently we have a self-consistent 
solution close to Tb (again we note that physically corresponds to graviational waves and will likely not affect 
the the dynamics close to the singularity). Suffice it to say that the numerical results discussed here indicate a 
self-consistent solution and serve to justify the dynamical arguments in [7]. 

In the case 7 < 4/3 we find numerically that /3" 1 — > 0, so that there is always an initial singularity as t — > — 00. In 
addition, we find that E1 1 , r, N x , £_ E x , all vanish as the initial singularity is approached (see FIGs. 6-8). 
However, the initial singularity is not, in general, isotropic. We shall discuss this, and the degenerate case 7 = 4/3, 
further in sections IV and V. 



IV. DERIVATION OF ASYMPTOTIC DYNAMICS 

From the numerical analysis we have the following conditions when 7 > 4/3: 

d: lim (E 1 1 ,B-\r,N^,N x ,^^ x ,^ + ,n,v,n u ,Q u ,n b -l) = 0, 

t— > — OO 

C 2 : d x (E 1 1 ,B- 1 ,r,N-,Nx,Y l -,Y.x,Y. + ,Q,,v,Q. u ,Q u ,Q,i > - 1) are bounded as -00. 
C3 : V — 0(f(t)) implies d x V = 0{f(t)) (asymptotic expansions in time 
can be differentiated with respect to the spatial coordinates). 

In particular, since Ex 1 — > as t — > —00, we can follow the analysis in [21] and use equations (107)-(124) to obtain 
the following asymptotic decay rates. 

Stage 1: First, C\ and C2 imply that q — ► 37 — 1. Using C\, C2 and the evolution equations, we obtain from 
Proposition 1 in [21] in succession, 



Ei 1 = 0( e K 3 T- 1 >- e ]*) 


(125) 


B- 1 = 0(e [37 - e]t ) 


(126) 


N x = O^t- 1 )-*]') 


(127) 


AT_ = 0(e^- 1 '>'^) 


(128) 


S x = ©(e 13 ^- 1 ^ 1 *) 


(129) 


£_ = Oic^- 1 ^) 


(130) 
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for any e > 0. It then follows from C\, C2, equation (126) and the algebraic expression for Of, that 

n = 0{c^- e ^) . (131) 

Stage 2: Using C\, C2 and C3, we obtain from equation (72) that 

r = C( c [( 3 T- 1 )- e ] t ) . (132) 

Using Ci, C2, C3 and the evolution equations, we obtain from Proposition 1 in succession, 

v = 0(e^-^-*) (133) 

£+ = ©(ePCT- 1 )-*]* + e^- 4 )^]') (134) 

Q u = o( e [( 3 ^ 4 )- e l*) (135) 

tt u = ^(el 2 ^- 4 )- 6 ] 4 ) . (136) 

Stage 3: First, 

g = 37-l + 5 (137) 

where g = 0(e^-^- e ^ + e^ 37-4 )- 6 !'), the dominant terms being £+ and v 2 . As in [21], we use Ci, C 2 , C 3 and 
Proposition 4 on the evolution equations to obtain 

El i =e^-V t [E 1 1 +0{g)} (138) 

B- 1 = e^B- 1 + O(g)} (139) 

iV x =c ( 3 7-i)*[7V x (140) 

N- = e^-V'lN- + O(g)} (141) 

S x =e 3 ^ t [t x +0(g)] (142) 

£_ = e 3(7 ~ 1)t [£_+0(.g)] (143) 

ft = c 37t [\/l2B- 1 +e>(; ? )], (144) 

where the hat variables are functions of x only. Then (72) gives 

r = e (*r-i)t[_B 1 i^ + 0(fl) ] . (145) 

The evolution equation for t> gives 

« = £e (37 ~ 4)t + ft , (146) 
where h = ©(e^ 37 " 1 )* + el 3 ' 37 " 4 ) -6 !'), the dominant terms being r and v 3 . The evolution equation for S + gives 

rS +e 3(7-l)*__|X_ 2 e 2(3 7 -4) t + OW | <7<2j7 ^| 

^+ I ^2 te 2t + S + e 2t + 0(«/i) 7=|- 
The evolution equation of Q u gives 

Q u = -27«e( 37 - 4 )' + Q u e 2 ^- 2 ^ + O(h) . (148) 
The evolution equation of Q u gives 

n u = -l^&J&i-W + Cl u c 2 ^- 2 ^ + 0(vh) . (149) 



We then eliminate the epsilons from g and h: 



0(e 3 (7-i)t + e 2 ( 37-4 )') |< 7 <2,7^ 
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g = 0(S + ) = ^ ^ 2t) - ^ 3 ^ '5- - ' ^ 3 (150) 
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and 

h = 0(r + v 3 ) = 0(c^-^ + c 3 ^- 4 )*) . (151) 
Summary: we have isotropization as t — > — oo for 7 > 4/3 with the asymptotic decay rates: 

E! 1 =e< 37 - 1 > t [£i 1 + 0(s)] (152) 

s -i =e 3 T t [B- 1 (153) 

r = eOT-iJtf.^i^ + . (154) 

AT X =e^- 1 ) t [^x+O(0)] (155) 

7V_ =e (37 - 1)4 [iV_+0( 5 )] (156) 

£ x =e 3 ^ 1 ) t [£ x +0( 5 )] (157) 

E _ = e 3 (T- 1 )*[S_+0( 5 )] (158) 

rE +e 3(7-l)*__|X_^ e 2(3 7 -4) t + OW | <7<2)7 ^5 

n = c 37t [a/125- 1 + C(. 9 )] (160) 

„ = ^ 37 - 4 )* + 0(h) (161) 

Q„ = -2 7 t)e^- 4 >< + Q u e 2 ^- 2 ^ + O(h) (162) 

Q u = -^W^-V* + n u e 2 ^- 2 » + 0(vh) , (163) 



(164) 



where 

a-OfS 1- /©(e^-^+e 2 ^- 4 )*) |< 7 <2, 7 ^f 
9 - ^ + >-\0(te 2t ) 7=| 

/i = 0(r + w 3 ) = ©(c^- 1 )* + c 3 ^- 4 '*) . (165) 

The asymptotic decay rates can be calculated numerically by plotting the log ratio of variables at different times. 
The case 7 = 1.8 is presented in FIGs 4-5. Note that the numerical calculations are consistent with the decay rates 
derived above. 

Next let us discuss the asymptotic dynamics when 7 = 4/3. Suppose the conditions C\, C2 and C3 again hold. 
Then, since Ei 1 — > as t — > —00, we can again follow the analysis in [21] and use equations (107)-(124) to obtain the 
asymptotic decay rates. Again, C\ and C2 imply that q — > 3, whence we obtain from Proposition 1 in succession, 

E 1 1 = 0(e( 3 - e )') (166) 

B- 1 = G(e (4 - e)t ) (167) 

N x = 0{c^- e)t ) (168) 

N- = 0(e (3 " e)t ) (169) 

S x = O^ 1 -^) (170) 

S_ = 0(e (1 - e ^) . (171) 

It then follows from C\, C2, C3, equations (72) and (167) that 

n = o(c( 4 - e )*) (172) 

r = 0(e (3 - e)t ) . (173) 

Unfortunately, v does not decay exponentially and we cannot continue with the integration method of [21]. 

However, we can still find the remaining decay rates heuristically. Equations (115) and (109) become asymptotically 

d t v = 2Y.+V + ... (174) 

9 t S+ = S+-^ 2 + ... , (175) 
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which give the approximate solutions 



(176) 



£+ = V + S+c( 1 - £ ) t + ... . (177) 

It follows that 

Qu = -l i v + Q u e {A - t)t + - (178) 
fl u = ~v 2 + fl u e^ t + ... . (179) 
These expressions are supported by the numerics. 

V. THE INVARIANT SET Si 1 = 

From numerical experiment we find that as t — ► — oo, Ei 1 — » very rapidly, 0, E_, E x , 7V X , JV_, r — > and 
_B — > co. We are now interested in the dynamics in the invariant set E\ l = at early times. In order to find the early 
time behaviour of fib (which is given in terms of fl, B and v by equation (96)) in the invariant set Ei 1 = 0, we need 
to obtain the evolution equation of fib from the evolution equations of ft, B and v as follows: 

dtflb = d v fl b d t v + dnflb d t fl + d B fl b d t B (180) 
where d t B, d t fl and d t v are given separately by (107), (114) and (115) with Ei 1 = 0. 

Equations (107)-(117) imply that Ei 1 = £_ = E x = N x = N- = ft = is an invariant subset of the invariant set 
Ei 1 = 0. In this invariant subset 

-\{Q b + Qu) = r = -^^-=Q (181) 

The dynamics in the invariant subset Ex 1 = £_ = E x = N x = N- = ft = (and hence r = Qb + Q u = 0) is given by 
J 2 + 6Q 6 + 6Q„ - 3 7 (1 + n 6 + Q„) 1 2 4 

+2 O b {(< + 1) - 3(1 - E +) (1 + 3(7 - 1)( ^; 2)+ ^ 2 ) + E+^j (182) 

dtv = - < 1 g V ^ ((2 - 7 )(l - 3E+) - 2( 7 - 1)) (183) 

1 ^ 3(2 7 -l)(l-r; 2 ) + 2 7 t; 2 ^ 1 

n (i + (2 7 -i)- 2 ) J" 6 + 3* 

+3E+ i 3 (1 + (27 - l)v 2 ) ^ + WT^-f u f 2 1 (l + (2 7 -l), 2 )G-G + | 7t 6 

-2 *{(,+ !)- 3(1 - E +) (1 + 3(7-D(l-^) + 7i> a ) + s+ |^ j (184) 

E+ = ^ (1 - fl b - ft u ) (185) 

1 , 1 ( 3(2 7 -l)(l-^ 2 ) + 2 7 .; 2 4 7 ^ 2 8^ 2 1 

9 -2 + 2\ (l + (2 7 -l), 2 ) " b + + (1 + (2 7 - l), 2 ) " b + (3T^) °"/ (186) 

Note that (187) 

Q u = -Q b = - TTT |^ TR O b (188) 



d t fi u = 2( q + 1) (fi b + ft u ) - 3(1 -v + )\fib + fi u + -{ y ^ /; o _ 1 ) n 6 + -a 



with 
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Since fib ~ fl 2 B 2 (see equation (96)), in this invariant set we have a closed system of ODE ((182)-(184)) in terms of 
the three variables fib, v and fl u . 

To obtain the equilibrium points, we let d t flb = 0, d t v = and d t fl u in equations (182)-(184). The results are 
summarized in TABLE I, where 



TABLE I. Equilibrium points and their stability 



equilibrium point 


eigenvalues 


stability 


6 = 0, v = 0, Q. u = 


-3(7-1), §7-L 1 


Saddle 


n b = o, v = o, n u = i 


-2(37- 2), 37-4, -1 


Saddle when 7 > 4/3. Sink when 7 < 4/3 


n b = i, v = o, n u = o (^i) 


37-4, 3(7 - 1), 2(37-2) 


Source when 7 > 4/3. Saddle when 7 < 4/3 


n b = o, u = ±i, n„ = o 


2, 4^i 


Saddle 


Q b = Fi( 7 ), v = ±F a ( 7 ), O u = F 3 ( 7 ) (P±) 


Ai( 7 ) > 0, A 2 (7) > 0, A 3 (7) > 


7 < 4/3. Source. 



TABLE II. Nonhyperbolic equilibrium point Tb = V± when 7 = 4/3 



The values of the variables 


Q b = 1, v = 0, n u = 0, Q b = Q u = E+ = 0, g = 3 


The eigenvalues and eigenvectors 


Ai = 1, vi = (1,0,0). A 2 = 4, = (-1,0,1). A 3 = 0, v 3 = (0, 1, 0) 
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p _ 1 j (37-4)(3 7 3 - 21 7 2 + 367- 16) _ / (3 7 -4)( 7 -l) (3 7 - 4)( 7 - 1)(9 7 2 - 28 7 + 16) 

1_ 3 7 (3 7 -2)(2- 7 ) 2 ' 2 "Y 2 7 2 - 7 7 + 4 ' 3_ 3 7 (3 7 - 2)(2 - 7 ) 2 

(189) 

For V±, we also obtain the values of q, £+, Qf, and Q„ from equations (185)-(188): 

« = i£t (190) 

^ + = ^-^(7) (191) 



n - + 2 (3 7 2 - 12 7 + 8) V (2 7 2 - 7 7 + 4) (3 7 - 4) ( 7 - 1) _ +p , , , . 

0b - ± 3 7 (37- 2)(2 - 7 ) 2 = 5(7) ( } 



n 2(37 2 - 127 + 8V(2 7 2 - 77 + 4)(3 7 -4) (7 - 1) 

0u = T 3 7 (37- 2)(2- 7 ) 2 =T 5(7) ( } 

Let us comment on the equilibrium points V±. -^2(7) is a real number only when 7 < 4/3, whence < f2j, < 1, 
> VL U > —1, \v\ < 1. So V± exist only when 7 < 4/3. The expressions for \, i = 1,2,3 are very complicated; to 
determine stability we calculated A^, i = 1,2,3 for different values of 7 < 4/3 numerically and found that Re(Xi) > 0, 
for i = 1, 2, 3 when 1 < 7 < 4/3. We conclude that V± are sources for the dynamics in the invariant subset given by 
Ei 1 = £_ = S x = N x = N_ = n = (= r = Q b + Q u ) when 7 < 4/3. 

Therefore in the invariant subset given by Ex 1 = £_ = £ x = N x = N- = ft = 0, T b is the global source for 
7 > 4/3, and V± (± depending on sign of v) are anisotropic sources for 7 < 4/3. There is a bifurcation at 7 = 4/3. 
Let us next discuss this degenerate case. When 7 = 4/3, V± and T b coincide. The values of the variables correspond- 
ing to this nonhyperbolic equilibrium point, the eigenvalues and the corresponding eigenvectors are given in TABLE 
II. At this nonhyperbolic equilibrium point, there exists a 2-dimensional unstable manifold and a 1-dimensional 
center manifold. The 1-dimensional center manifold is tangent to t>3 = (0,1,0) at this point, therefore the center 
manifold has the form tt b - 1 = a Q v 2 + a x v 3 + a 2 v 4 + a 3 v 5 + 0(v 6 ) and tt u = b Q v 2 + b x v 3 + b 2 v 4 + b 3 v 5 + 0(v 6 ). 
The stability is determined by the dynamics in the center manifold. We find that the center manifold is g iven by 
il b - 1 = -ft; 2 - ±f-v 4 + 0{v 5 ) and fl u = -§v 2 + ^-v 4 + 0(v 5 ) and that the dynamics in the center manifold is 
consequently determined by ^ = yu 3 + ^jf-v 5 + 0(v e ). From the dynamics in the center manifold, we find that the 
nonhyperbolic equilibrium point is a source for the center manifold, therefore the nonhyperbolic equilibrium point is 
a source for the 3-dimensional invariant subset Ei 1 = X_ = S x = iV x = iV_ =O = 0. This is consistent with the 
decay rates calculated above and with numerical analysis (and phase portraits in the 3-dimensional invariant set). 
Therefore models isotropize 'slowly' in the radiation case. In summary, when 7 > 4/3, T b is a global source; when 
7 = 4/3, T b = V± is still a global source; when 7 < 4/3, T b becomes a saddle and a new pair of equilibrium points 
V± appear as a pair of sources. In addition to Tb and P±, there are also a number of saddle equilibrium points (see 
TABLE I). 

We claim that when 7 > 4/3, T b (Ei 1 = £_ = S x = N x = iV_ = fl = r = Q b = Q u = = v = 0, fl b = 1) is also 
a source for the full state space and that when 7 < 4/3, V± (Ei 1 = £_ = S x = N x = iV_ = fl = r = S + = 0, O5 = 
^1(7), v = ±F 2 (7), Q u = F 3 ( 7 ), S+ = F 4 ( 7 ), Qb — ±^5(7), Q u — TFsil)) are sources for the full state space. We 
give an argument for our claim as follows. First, we observe that at Tb 

q + 1 = 3 7 > 0, q + 3S+ - 2 = 3(7 - 1) > 0, Ti = I > 0, T 2 = 3 7 - 1 > 0. 



and at V± 



27 2(7-1) (37-2)7 2 

q+l = —^>0, g + 3S + -2 = ^-^ >o, r 1 = - V 11 >0, T 2 = - >0. 

2 — 7 (2 — 7) (27 2 -77 + 4) 2 — 7 



where 



Ti = (1 + (7 ^ ^ {q + 1) - \{l - 3E+)(1 + v 2 ) - 1, T 2 = 2(q + 3S+ - 1) - |(-n 6 - 0„ + P b + p u + n fe + n„). 
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So near the equilibrium points Tb and V±, q + 1 > and q + 3X + — 2 > and thus as t — > — oo, we have B — > +oo and 
Pi 1 — > 0. Near JF b and P±, we can neglect the terms with l E\-d x ' in the equations (107)-(116) and treat the PDE as 
a system of ODE (see section IV). Then from the linearization of equations (110)-(114) near Tb and V± and using 
the fact that q + 3S+ — 2 > 0, T\ > near Tb and V±, we find that as t — > — oo, S x , iV x , AL. and (therefore 
Q) will decrease monotonically towards zero. We also observe that near Tb and V± 

d t (Qb + Qu) = 2(q + 3S+ - l)(Q b + Q u ) 

+ ^3(at_e x - v x s_) - ^(Q + Q b + Qu?j (-n b -n u + p b + Pu + iib + n u ) (194) 
= (2( q + 3S+ - 1) - |(-n 6 _ q u + p b + p u + n 6 + n„)^ (q 6 + q u ) 
+ ^3(iv_s x - tv x s_) - ^ (_n 6 - + p fc + p„ + n fc + n„) . (195) 

Linearizing this evolution equation for Qf, + Q u near Tb and P± and using T 2 > 0, we conclude that Qb + Qu, and 
thus r (= 3(-/V_£ x — ]V X E_) — §(Q + Q/, + Q u )), decrease towards zero as t — > — oo. 

In summary, near P{, and P±, Pi 1 , E_, S x , AT X , N-, and r will decrease monotonically towards zero as t — > — oo, 
so all orbits in the full state space evolve back towards the invariant subset given by Pi 1 = £_ = S X = 7V X = 7V_ = 
O = 0. Because Tb (when 7 > 4/3) and V± (when 7 < 4/3) are sources in this invariant subset, we expect that the 
orbits mentioned above will shadow orbits in the invariant subset and evolve back towards Tb (when 7 > 4/3) or V± 
(when 7 < 4/3). Therefore, Tb is a source for the full state space when 7 > 4/3 and V± are sources for the full state 
space when 7 < 4/3. The equilibrium points V± correspond to new anisotropic brane world cosmologies. Note that 
the tilt is not extreme (v 2 ^ 1) at V±. 

The numerical experiments have confirmed that Tb is a source in the full state space when 7 > 4/3 (see FIGs.1-3). 
The numerical simulation is consistent with the decay rates given in section IV and with the analysis given above 
that Tb is a source for the full state space when 7 = 4/3. Therefore models isotropize 'slowly' in the radiation 
case. Next, we present numerical evidence that the V± are sources when 7 < 4/3. For example, when 7 = 1.2, 
il b = 0.767361111, v = ±0.2294157339, O u = -0.1006944444 (and Ai, 2 = 0.27236 ± 0.4838767293i, A 3 = 2.644927); 
in this case q = 1.9999999, £+ = 0.1666666, Q h = ±0.3935117 and Q u = T0.3935117. These results for 7 = 1.2 arc 
consistent with the results given by numerical experiment (see FIGs.6-8). 



VI. DISCUSSION 

All models have an initial singularity as t — > —00. In addition, we find that {O, £_, iV x , S x , N_, r} — > as t — > —00 
for all 7 > 1. In the case 7 > 4/3, the dynamical and numerical analysis indicates that {v, Sl u , Q u , £+} — > (and 
il b — » 1) for all initial conditions. In the case of radiation (7 = 4/3), the models still isotropize as t — > —00, albeit 
slowly. For 7 < 4/3, {v, fl u , Q u , E + } tend to constant but non-zero values as t — > —00. 



A. Tilt 

In the invariant set v = 0, all models isotropize to the past (for 7 > 1). Thus in the spatially homogeneous case 
with no tilt (with Pi 1 = 0), it follows that there exists an isotropic singularity in all orthogonal Bianchi brane- world 
models in which Of, dominates as the initial singularity is approached into the past (consistent with the results of 
[11]). In particular, T b is a local source and in general the initial singularity is isotropic. The linearized solution 
representing a general solution in the neighbourhood of the initial singularity in a class of Gi models was given in [7] ; 
it was found that Tb is a local source or past-attractor in this family of spatially inhomogcncous cosmological models 
for 7 > 1. The exponential decay rates of the G2 models (about T b ) given in [7] are consistent with those here in 
v = case (for all 7). 

Let us now assume that 11 / (the general case). For the Gi models in the timelike area gauge we recover the 
orthogonally transitive tilting Bianchi type VIo and VIIo models in the spatially homogeneous limit (with one tilt 
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variable). (A spatially homogeneous cosmology is said to be tilted if the fluid velocity vector is not orthogonal to the 
group orbits, otherwise the model is said to be non-tilted). There are no sources in the Bianchi type Vlo and Vllo 
models with tilted perfect fluid; the past attractor is an infinite sequence of orbits between Kasner points [22]. A 
description of the dynamics of tilted spatially homogeneous cosmologies of Bianchi type II (with a perfect fluid with 
linear equation of state) has been presented [23]. The Bianchi II cosmologies, while very special within the whole 
Bianchi class, play a central role since the Bianchi II state space is part of the boundary of the state space for all 
higher Bianchi types (including types Vlo and VIIo). The class of tilted Bianchi II cosmologies can be described by 
a set of expansion-normalized variables, and the state space is bounded. (Note that expansion-normalized variables 
and (3 -normalized variables are effectively the same close to the initial singularity). In more detail [23], there is no 
equilibrium point that is a local source, except in the special case 7 = 2 (in which there is a local source, namely a 
subset of the Jacobs disc). In the tilting situation there are two Kasner circles, the standard Kasner circle and the 
Kasner circle with extreme tilt. The evolution of tilted cosmologies of Bianchi type II in the singular asymptotic 
regime is governed by infinite hcteroclinic sequences which contain orbits that join two points between these two sets. 
This is analogous to the case of non-tilted spatially homogeneous cosmologies of Bianchi types VIII and IX which 
exhibit so-called Mixmaster oscillatory behaviour as the singularity is approached into the past. Earlier work had 
shown that there are no local sources in Bianchi type Vlo models with a magnetic field, so that in general these 
models also exhibit mixmaster behaviour to the past [24]. 

Consequently, the Bianchi models have bifurcations at 7 = 2/3 and 4/3 (e.g., the dimension of the unstable and 
stable manifolds of the equilibrium points change). There is no local source in general, and the models exhibit 
mixmaster behaviour to the past. A subset of models are past asymptotic to the flat Friedmann model (i.e., the 
singularity is isotropic) if 7 > 4/3. In particular, there is a bifurcation at 7 = | in both the tilting and magnetic field 
models in general relativity. We note that this is consistent with our bifurcation 7 = 4/3 (with 7 — > 27 in brane- world 
models). 

B. Extensions 

The earliest investigations of the initial singularity, which used only isotropic fluids as a source of matter, suggested a 
matter-dominated isotropic singularity for all 7 > 1 [6,11,7,25]. However, it was shown in later work using anisotropic 
stresses [9] that the initial singularity for a magnetic brane-world could be cither locally isotropic or anisotropic. In 
particular, Barrow and Hervik [9] studied a class of Bianchi type I brane-world models with a pure magnetic field 
and a perfect fluid with a linear barotropic 7-law equation of state. They found that when 7 > |, the equilibrium 
point Fb is again a local source (past-attractor), but that there exists a second equilibrium point denoted PH\, which 
corresponds to a new brane-world solution with a non-trivial magnetic field, which is also a local source. When 7 < |, 
PH\ is the only local source. This was generalized by [10], who presented a thorough investigation of the initial 
singularity in brane-world cosmological models. It was shown that for a class of spatially homogeneous brane- worlds 
with anisotropic stresses, both local and nonlocal, the brane-worlds could have either an isotropic singularity or an 
anisotropic singularity; indeed, using a continuity argument it was shown that there exists a past attractor for models 
with nonlocal anisotropic stresses of type Vp. v = UD^ U where d t D^ u is sufficiently small. Hence, there is a class of 
models with V^ v 7^ which have an anisotropic past attractor. How large this class is, and if this anisotropic past 
attractor exists for generic brane-worlds, needs further work. However, the analysis in the present paper is consistent 
with the results of [9,10] in which an isotropic singularity exists for 7 > 4/3. 

In general, is not specified. It must be derived from the exact 5-dimensional field equations in a self-consistent 
way. In our analysis we have assumed that the effective nonlocal anisotropic stress is zero (in the fluid comoving 
frame). Indeed, this is the only assumption we have made. But it is expected that inclusion of will not affect 
the qualitative dynamical fatures of the models close to the initial singularity (and we still expect isotropization at 
early times). £ M „ can be irreducibly decomposed according to equation (6). Hence, we expect that V^v ~ Ug^v on 
dimensional grounds, and so for a Friedmann brane close to the initial singularity we expect that V^v ~ a 2 UC^ 
(where is slowly varying) , which is consistent with the linear (gravitational) perturbation analysis (in a pure AdS 
bulk background) [26]. Hence, V^v is negligible dynamically close to the initial singularity. 

The results might also be applicable in a number of more general situations. For example, in theories with field 
equations with higher-order curvature corrections (e.g., the four-dimensional brane world in the case of a Gauss- 
Bonnet term in the bulk spacetime [27]), the results concerning stability are not expected to be affected since the 
curvature is negligible close to the initial singularity. 
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VII. CONCLUSIONS 



Therefore, the numerical analysis supports the fact that in spatially inhomogcncous G2 brane-world cosmological 
models the initial singularity is isotropic [8]. 

Therefore, unlike the situation in general relativity, it is plausible that typically the initial singularity is isotropic 
in brane world cosmology. Such a 'quiescent' cosmology [28], in which the universe began in a highly regular state 
but subsequently evolved towards irregularity, might offer an explanation of why our Universe might have began its 
evolution in such a smooth manner and may provide a realisation of Penrose's ideas on gravitational entropy and 
the second law of thermodynamics in cosmology [29]. More importantly, it is therefore possible that a quiescent 
cosmological period occuring in brane cosmology provides a physical scenario in which the universe starts off smooth 
and that naturally gives rise to the conditions for inflation to subsequently take place. 

Cosmological observations indicate that we live in a Universe which is remarkably uniform on very large scales. 
However, the spatial homogeneity and isotropy of the Universe is difficult to explain within the standard general 
relativistic framework since, in the presence of matter, the class of solutions to the Einstein equations which evolve 
towards a Robertson- Walker universe is essentially a set of measure zero. In the inflationary scenario, we live in an 
isotropic region of a potentially highly irregular universe as the result of an expansion phase in the early universe 
thereby solving many of the problems of cosmology. Thus this scenario can successfully generate a homogeneous and 
isotropic Robertson- Walker-like universe from initial conditions which, in the absence of inflation, would have resulted 
in a universe far removed from the one we live in today. However, still only a restricted set of initial data will lead to 
smooth enough conditions for the onset of inflation. 

Let us discuss this in a little more detail. Although inflation gives a natural solution of the horizon problem of 
the big-bang universe, inflation requires homogeneous initial conditions over the super- horizon scale, i.e., it itself 
requires certain improbable initial conditions. When inflation begins to act, the universe must already be smooth 
on a scale of at least 10 5 times the Planck scale. Therefore, we cannot say that it is a solution of the horizon 
problem, though it reduces the problem by many orders of magnitude. Many people have investigated how initial 
inhomogeneity affects the onset of inflation [30,31]. Goldwirth and Piran [30], who solved the full Einstein equations 
for a spherically symmetric spacctime, found that small-field inflation models of the type of new inflation is so sensitive 
to initial inhomogeneity that it requires homogeneity over a region of several horizon sizes. Large-field inflation models 
such as chaotic inflation is not so affected by initial inhomogeneity but requires a sufficiently high average value of 
the scalar field over a region of several horizon sizes [32]. Therefore, including spatial inhomogeneities accentuates 
the difference between models like new inflation and those like chaotic inflation; inhomogeneities further reduce the 
measure of initial conditions yielding new inflation, whereas the inhomogeneities have sufficient time to redshift in 
chaotic inflation, letting the zero mode of the field eventually drive successful inflation. In conclusion, although 
inflation is a possible causal mechanism for homogenization and isotropization, there is a fundemental problem in 
that the initial conditions must be sufficiently smooth in order for inflation to subsequently take place [11]. We have 
found that an isotropic singularity in brane world cosmology might provide for the necessary sufficiently smooth initial 
conditions to remedy this problem. 

It would be of interest to study general inhomogeneous (Go) brane world models. The exponential decay rates in 
the case 7 > 4/3 are calculated in the Appendix (in the separable volume gauge using Hubble- normalized equations). 
The decay rates are essentially the same as in the G2 case studied in section IV (there are minor differences due to 
/3-normalization and the absense of two tilt variables in the G2 case: cf. equations (157)-(159)). This supports the 
possibility that in general brane world cosmologies have an isotropic singularity. We hope to further study Go brane 
world models numerically in the future. 
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APPENDIX A: HUBBLE-NORMALIZED G EQUATIONS 

We introduce _ff-normalizcd variables as follows: 
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(AT 1 , Ej) := (N~\ eJ)/H 
(U a , A a , £ Q/3 , N a(3 , R a ) := (u„, a", a Q/3 , n a/3 , 

6k 2 6k 2 
(O, fi u , Q u ) := (k 2 p, — p„, — q u )/(3H 2 ) 



(Al) 
(A2) 

(A3) 
(A4) 



Using the separable volume gauge: U a = r a , with Af = 1, as the choice of temporal gauge, and the Fermi-propagated 
gauge: R a = 0, as the choice of spatial gauge, the equations are: 



where 



a/3 



d t Ej = (qS a 
d t H = -{q + l)H 
d t r a = {q Sj 3 - Y.J 3 ) rp + d a q 

d t A a = (qS a -J: a (} )A' 3 + ^d^ 

d t Z af3 = (q - 2) S a/3 - 2JV <a 7 N^i + N {a(3) - <F< Q (0 7 - r 7 ) A# 

+ (a 7 - 2A 7 ) A^a + (<F< Q 9 7 + A^ a ) + 311^ 
d t N al3 = {q 6 {a s + 2£ (a 5 ) N^ 5 - e' 5{a 9 7 £« 4 

d t n = -^-v a d a n + g; 1 [ 2G +<Z - (3 7 - 2) - (2 - 7 ) v 2 - 7 (£ a/3t /V) 

- 7(9 a - 2A Q )w a +71;" d Q lnG+] ft 

<9 t u Q = - 1/ 3 9/31;" + S af3 dp In G+ - 12L_12 (i _ ^ ^/s ^ Jn n _ 2 ^ 

7 

+ GI 1 [( 7 -l)(l-« 2 )(a (3 ^)-(2- 7 ) W ' 3 9 (3 lnG + 



+ 



(7-1) 



(2 - 7) (1 - w 2 ) v* 3 {dp InO - 2^) + (3 7 - 4) (1 - w 2 ) 



+ (2 - 7) (£ /37 «'V) + G_ (r^) + [G+ - 2(7 - 1)] (Apv 13 ) 
- T. a p v 13 - r a - v 2 A a + e a/37 N 0S u 7 v s 

d t n u = 2( q - i)Q, 



-—^ af3 v a v p n u + 2A a Q a u - d a QZ + v a X a 

O "T" V 



d t Q a u = 2(q - 1)QZ - ^ a pQ f3 + 4 — ^ — A n u 

+ -^[2v a v? - (1 + v 2 )S^}r fj Q u - dp{P u 5^ + Ilf) + X c 



x a = 



j 2 (l~v 2 ) 



6G+ 3 G 
+ y 2 (l-v 2 ) 



M af3 d p v 2 



2^2 



7(1 -v 2 ) 
6G+ 2 G_ 



■M af3 n 2 H 2 [d lnft-2773] 



[o^ + (3 - w 2 ) - x^v" - 2A^]v a n 2 n 2 



6G+ 2 G 

M Q/3 = G-5 a(3 + (7 - l)^^ 

9 = 2S 2 + l(O tot + 3P tot ) - i (0 Q - 2A a ) r° 
r a = —d a \nH . 



The constraints are: 



= 2 (8 [a - r [a - A [a ) Erf - e af3S N 5 ~< Ej 
= 1 — flk — S 2 — O tot 

= dpir* 1 + {2S a f3 - Y, a p) r 13 - ZA fj S Q/3 - e a ^ N fjS £,/ + 3Q c t 



(A5) 
(A6) 
(A7) 

(A8) 



(A9) 
(A10) 

(All) 



(A12) 
(A13) 

(A14) 



(A15) 

(A16) 
(A17) 
(A18) 

(A19) 
(A20) 
(A21) 



21 



where 

ftfc = - I (2d a - 2r a - 3A a ) A a + l - (N a0 N a ^) - -L (N a a ) 2 . (A22) 

If we have (for 7 > 4/3): 

C{: lim (E a \H-\r a ,A a ,Z a P,N af:) ,n,v a ,n u ,QZ,n b -l) = , 

t— > — CO 

C* 2 : di{E a \U~ x ,r a , £ a/3 , 7V a/3 , ft, ft u , ft fc - 1) are bounded as i -» -00. 
C3 : V = 0(f(t)) implies 9jF = 0(f(t)) (asymptotic expansions in time 
can be differentiated with respect to the spatial coordinates). 

then we can follow the analysis in [21] to obtain the asymptotic decay rates. 

Stage 1: First, C\ and C\ imply that q — > 37 — 1. Using C*, and the evolution equations, we obtain from 
Proposition 1 in [21] in succession, 





(A23) 




(A24) 


r a = 0(e^-V~*) 


(A25) 


A a = 0(e [(37_1) - e]t ) 


(A26) 


N al3 = C(e [(37_1)_elt ) . 


(A27) 



It then follows from C^, C£, equation (A24) and the algebraic expression for ft;, that 

ft = 0{e [ ^-~ e]t ) . (A28) 
Stage 2: Using C{, Cj, C3 and the evolution equations, we obtain from Proposition 1 in succession, 

S"/ 3 = 0( e [3(7-l)-e]* + e [2(3 7 -4)-e]t^ 

Stage 3: First, 

q = 3 7 - 1 + <D(e^-V-*) . 
As in [21], we use C* , C£, C3 and Proposition 4 and the evolution equations to obtain 

Ej = e^-V t [E a i + 0(g)} 
H' 1 = e^H' 1 + 0{v 2 )] 
r a =e^- 1 »[r a + 0(g)] 

A a =e ( 3 T- 1 )*[i Q + 0( 5 )] 
=e 3(7-l)*[ i V«/3 + C ,( ff )] 

= e 3 T*[\/l2W- 1 +O(u 2 )] . 

s «/3 = J 3^5«< a ^>e 2 ( 3 T- 4 )* + E^e 3 ^" 1 )' + I < 7 < 2, 7 ^ f 

I lOfi^t^te 2 * + t a(3 e 2t + Q{vh) 7 = f 

QZ = -2 7 v a e^-^ + Q2e 2(37 " 2) * + 0(h) 

n u = -^ 7 w^- 4 )< + n u e 2 ^- 2 ^ + o(vh) 
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(A29) 
(A30) 
(A31) 
(A32) 

(A33) 

(A34) 
(A35) 
(A36) 
(A37) 
(A38) 
(A39) 
(A40) 

(A41) 

(A42) 
(A43) 



where g = 0{e^~^- f ^ + e [ 2 ( 3 T- 4 )- £ ]*), the dominant terms being E Q/3 and v 2 ; and h = ©(e 3 ^- 1 )* + e [ 3 ( 3 T- 4 )- e l*). 
Finally we eliminate the epsilons from g, h and other O terms by repeating (and using the results from) Stage 3: 

g = 0(E) = | ^ e2t) +C )3 7 1 7 | ' 77S (A44) 

/i = 0(r + v 3 ) - 0(e^-^ + e 3 ^" 4 )') (A45) 
g = 3 7 - 1 + 0(e 2 ^- 4 ) 4 ) (A46) 

Summary: as t — > — oo we have that 

= e^- 1 )*^* + 0(g)] (A47) 

W" 1 =e 37 *[W- 1 +0(w 2 )] (A48) 

r a =e^- 1 ) t [r a +0( 5 )] (A49) 

A a = e^ 7 " 1 ^! + 0{g)] (A50) 

N al3 = e 3 (T^ 1 )*[iV Q ' 3 + O(g)} (A51) 



where 



~ 1 Wv^vPhe 2 * + £ a Pe 2t + <D{vh) 7 =| 



(A52) 



n = e 37 *^^ -1 + C(w 2 )] (A53) 

/ = v a e^- i ^ t + 0{h) (A54) 

QZ = -2 7 t) Q e (37 - 4)t + Q2e 2(37_2) * + 0(h) (A55) 

Q„ = -^W^" 4 * + tt u e 2 ^- 2 ^ + 0(vh) , (A56) 



.-^>-{$£ w+ ^ )1 7 < -V J,T '' 1 (A57) 

h = 0(r + u 3 ) = 0(e (37 - 1)4 + e 3(37 - 4) *) . (A58) 
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FIG. 5. Exponential decay rates: 7 = 1.8, ti 
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FIG. 7. Initial singularity for 7 < 4/3: 7 = 1.2, £ = -40 
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FIG. 8. Initial singularity for 7 < 4/3: 7 = 1.2, £ = -90 
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